Collective vibrational states with fast iterative QRPA method 
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An iterative method we previously proposed to compute nuclear strength functions [lj] is developed 
to allow it to accurately calculate properties of individual nuclear states. The approach is based 
on the quasi-particle-random-phase approximation (QRPA) and uses an iterative non-hermitian 
Arnoldi diagonalization method where the QRPA matrix does not have to be explicitly calculated 
and stored. The method gives substantial advantages over conventional QRPA calculations with 
regards to the computational cost. The method is used to calculate excitation energies and decay 
rates of the lowest lying 2 + and 3~ states in Pb, Sn, Ni and Ca isotopes using three different Skyrme 
interactions and a separable gaussian pairing force. 

PACS numbers: 21.60.Jz, 21.10.Re 
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I. INTRODUCTION 

The goal of nuclear structure theory is to be able to 
predict and model the physics of the atomic nucleus. This 
involves the ground-state properties, as well as different 
modes of excitation and decay. One of the possible meth- 
ods to compute excited states in nuclei is based on the 
quasi-particle-random-phase approximation (QRPA) |2j. 
This approach can be derived by considering the linear 
response of a nucleus when perturbed by an external field. 
From the response one can extract information about ex- 
cited nuclear states and cross sections for nuclear reac- 
tions. The QRPA approach is particularly interesting in 
connection with nuclear density-functional theory (DFT) 
as the method can be applied also when starting from a 
density functional. In order for the QRPA method to be 
practical, it is very important to implement it in ways 
that have low computational costs. For phenomenolog- 
ical DFT approaches, a low computational cost would 
allow dynamical properties to be considered when fine 
tuning model parameters. A numerically efficient method 
is also essential for applications to deformed and heavy 
nuclei which are otherwise prohibited by the time and 
memory required to construct and diagonalize the large 
QRPA matrix. 

Two recent solution methods address these issues. The 
Finite Amplitude Method (FAM) d, 0] generates the re- 
sponse of a nucleus to an external field by solving the 
linear response equations iteratively for each requested 
external field frequency. FAM furthermore uses the same 
mean fields as in the Hartree-Fock-Bogoliubov (HFB) 
ground-state calculation and employs finite differences 
to linearize the equations of motion. In its current form 
FAM uses a smoothing method to improve stability and 
therefore one cannot easily extract the exact QRPA eige- 
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namplitudes. The same is true for the iterative Arnoldi 
method [l| which is able to provide smoothened QRPA 
strength functions and their energy weighed moments, 
but does not generate accurate individual states. A com- 
mon aspect of both methods is however their ability to 
generate partial solutions of the full QRPA problem with 
reduced computational effort. 

The purpose of this paper is to generalize the Arnoldi 
method which we previously developed for iterative cal- 
culations of RPA strength functions. The generalization 
involves nrodifying the method so that it becomes possi- 
ble to not only compute strength functions but also sets 
of individual excited states with high accuracy. As a first 
step the new method is applied to the calculation of exci- 
tation energies and decay rates of the lowest lying 2 + and 
3~ states in several isotope and one isotone chain. Par- 
ticular focus is given to the region around double-magic 
208 Pb where new experiments are currently planned [5|. 

This paper is organized as follows: in Sec. II the QRPA 
formalism is briefly reviewed and specific aspects of our 
formulation are discussed. In Sec III the computational 
cost and accuracy of the method is evaluated. In Sec. 
IV the method is applied to the calculation of energies 
and transition probabilities of the lowest J n = 2 + and 
J n = 3~ states in a selection of semi-magic even-even 
nuclei. Finally conclusions are given in section V. 



II. QRPA IN TERMS OF FIELDS 

The iterative method is based on the QRPA equations 
[2J, lfj-l8| which can be derived by starting from time- 
dependent HFB theory. Here we present the main parts 
of the derivation, highlighting aspects relevant to the it- 
erative formulation. In cases where the expressions are 
not fully defined we use notation consistent with Ref . [2J . 

The QRPA equations can be derived by considering 
a general time-dependent wavefunction which is oscillat- 
ing between the ground state and an excited state with 
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We limit the consideration to small-amplitude oscilla- 
tions around the ground state so that the corresponding 
generalized density-matrix TZ [2| can be expanded to first 
order in C exc 



This expansion is taken to first order in the transitional 
fields, which is enough for small-amplitude vibrations and 
leads to 
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In order to make use of the time-dependent-Hartree- 
Fock-Bogoliubov (TDHFB) equations of motion, it is 
desired that the time-dependent density should be a 
HFB density at all times (i.e., it should be a projector 
TZ 2 = TV). We further assume that TZ gs can be approx- 
imated with the ground-state HFB density. Then the 
most general approximation for the transition density 
which ensures that TZ (t) a projector for small-amplitude 
vibrations involves both the forward Z and the backward 
Z'* going amplitudes 
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In this expression we have made use of the matrix 
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In our case with a density-independent pairing interac- 
tion it is only the h = d£ /dp field which is non-linear 
in the densities and becomes linearized. With a density- 
dependent pairing interaction the A, A' fields would also 
have to be linearized and would give an additional con- 
tribution to the h field. 

Inserting the expressions for fields and densities into 
the QRPA equation, and multiplying from the left with 
W and from the right with 14, gives a system of equations 
for the unknown excitation energies hu and the Z and Z' 
amplitudes: 
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written in terms of the U and V pairing matrices [2j re- 
lated to the HFB ground state as well as the matrices 
of quasiparticle operators [aa'J = ctiaj. It should be 
noted that if the wavefunctions in the beginning were 
taken as HFB vacuums one would not obtain any back- 
ward going amplitudes, as can be seen from Eq. ([T]) by in- 
serting the HFB ground state. Thus it is the assumption 
of the transition density TZ being as general as allowed 
by the TZ 2 = TZ criteria, which allows for the existence of 
the implicitly defined correlated ground state. 

Inserting the expression for TZ(t) into the TDHFB 
equations of motion ih -4 - = [H,TZ(t)] [8| and taking 
the small-amplitude limit leads to the QRPA equation 



fruTZ- [H[TZ gs },TZ] + [H^iZlTZgs]. 

In this expression the hermicity property Ti l \TZ] = 
T-L 1 \TV] ] of the effective interaction was assumed and 



the time-dependent fields Ti(t) 
around the ground-state value 
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hwZ = EZ + ZE + W, 

-TuuZ^ = EZ'^ + Z^E + W n . 

In this equation, E denotes a diagonal matrix of positive 
quasi-particle energies and the W matrices depend on the 
linearized fields 

W = U r hV* + U^hU* + V^h'^V* -V r h T U*, 
W ,Jf = V T hU + V T AV + U T A" i U -~U T h T V, (2) 

which can be expressed in terms of the transition densi- 
ties 

p = uzv T + v*z /J{ u\ 
k = UZU T + v*z^v\ 
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It is instructive to look back and consider the approx- 
imations used in the derivation of these equations. The 
main approximations appear to be the use of the TDHFB 
equations of motion, which restricts us to the consid- 
eration of excited states connected by two quasiparticle 
operators to the ground state and the assumption that 
the ground-state density can be approximated with the 
density of the HFB ground state. 



As an example to illustrate the iteration procedure, 
we neglect spin and isospin and consider a term in the 
energy of the form 

£[p]=Jp a+2 (r)dr, 

which gives the linearized field 

h im = (a + 2) (a + 1) J 4>l (r) {p a g . s . (f) P (f)) c/> m (r) dr. 

In this case, the action of the QRPA matrix on an eigen- 
vector can be calculated in three steps. 

The first step is to generate the densities p according to 
Eq. $5§ and expressing them in r space. For the next step, 
h is calculated as above. Alternatively, in the case of a 
density-independent interaction, where fields are already 
linear in densities T-L l \K] = H[1Z], this can be achieved 
using the HFB mean-field routines for calculating matrix 
elements. Finally multiplying the fields with U and V 
matrices as in Eq. $Z§ one obtains the W matrices. 

The main advantage of expressing the equations in this 
form is that calculating and storing two-body matrix ele- 
ments can be avoided and instead one can rely on the ex- 
pressions for HFB fields. The price to pay is that the den- 
sities and integrals for the matrix elements of the fields 
are recalculated for each matrix vector product, in the 
same way as when performing the HFB iterations to find 
the ground state. Thus it is important to investigate how 
many iterations i.e. matrix-vector products are needed 
in order to obtain acceptable convergence, and whether 
the iteration procedure introduces numerical errors which 
could lead to instabilities. 



III. ACCURACY AND EFFICIENCY OF THE 
METHOD 

The iterative QRPA solver is implemented by extend- 
ing the program hosphe (vl.02) [9| and will be included 
in the next published version of the program. This code 
uses a spherical harmonic oscillator basis and takes ad- 
vantage of the Wigner-Eckart theorem in order to work 
with angular momentum reduced quantities. The use of 
reduced quantities keeps the HFB and QRPA dimensions 
small and makes the code a useful tool for testing differ- 
ent calculational methods. 

In order to verify that the QRPA implementation is 
correct, a comparison is made with a recent QRPA im- 
plementation based on the HFBTHO code |1(J]. This code 
is able to treat axially deformed nuclei and its QRPA im- 
plementation is based on the traditional diagonalization 
of a large QRPA matrix. Therefore, applications of this 
code are limited to cases where dimensions can be kept 
within manageable limits. 

As a test case we consider the nucleus 1 |Oio and com- 
pare the ground-state energy and energies of the QRPA 
excitations obtained in both codes. In order to have a 
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Table I: Comparison of HFB and QRPA calculations per- 
formed for the nucleus ls O without any pairing truncation 
and without Coulomb interaction. For the Skyrme interaction 
we use the SLy4 parametrization [111 ] with a delta (volume) 
pairing interaction |T3 | with strength Vo = —200 MeV(fm) 3 
and a one-body center of mass correction. The results are 
obtained with a harmonic oscillator basis where the maxi- 
mum oscillator shell included has principal quantum number 
Nmax — 5 and the oscillator constant is set to 0.865 (fm)~ . 
In this table, the transition strengths are calculated using the 
isoscalar transition operators of Ref. (l4l | . Energies have units 
of MeV and the B (EI) transitions are in units of e 2 (fm) . 
For each multipolarity the state lowest in energy with an ap- 
preciable strength is compared. When both codes give the 
same decimals, they are printed in bold. 



benchmark result that is useful for testing future QRPA 
codes we list values obtained from both codes in Table [I] 
Several different recipes on how to truncate the pairing 
space and how to treat the Coulomb interaction exist in 
the literature, so in order to make the benchmark results 
as useful as possible, the results are obtained without any 
pairing truncation and without any Coulomb interaction. 
The remaining parameters are listed in the caption of 
Tab.U 

The implementation based on the axially deformed HF- 
BTHO code [H| allows us to test its accuracy by perform- 
ing calculations for excitations with different angular mo- 
mentum projections on a principal axis of the nucleus. 
Since the comparison is made for a spherical nucleus, 
these different calculations should ideally give the same 
result. In this way, for the QRPA implementation of 
[10| . the precision of the 2 + excitation was estimated to 
be roughly 10~ 4 both for the energy and for the B (E2) 
value. As seen from Table []] the lowest states calculated 
with both codes agree to about this precision. A similar 
accuracy test with the HOSPHE code is not possible as 
it works in spherical symmetry. But since HOSPHE uses 
the same mean-fields both in the QRPA and the HFB 
calculations we expect about the same accuracy for the 
QRPA excitations as for the groundstate energy. The full 
strength functions calculated with both codes were also 
compared and turned out to be indistinguishable by eye 
when the reduced transition probabilities are plotted as 
function of the energy of the excited states. 



A. Iterative solutions 

As described above, the product of the QRPA matrix 
acting on an arbitrary vector can be calculated with- 
out constructing the matrix explicitly. When this tech- 
nique is used, traditional matrix diagonalisation rou- 
tines, which need explicit information about the ma- 
trix elements can not be used. Instead, one must re- 
sort to indirect iterative methods, such as the Lanczos 
or Arnoldi [17] methods. For non-hermitean problems, 
such as the QRPA eigenvalue problem, the Implicitly 
Restarted Arnoldi method (IRA) [18|, [jjj] is one of the 
most commonly used methods for finding accurate ap- 
proximations for the eigenstates lowest in energy. IRA is 
a more advanced version of the original Arnoldi method, 
giving faster convergence and a reduction in computa- 
tional effort. As with the iterative Arnoldi method of 
Pi, the IRA method generates a set of basis vectors, usu- 
ally called Ritz vectors, which span a vector space called 
the Krylov subspace, and uses these vectors to represent 
the QRPA eigenvectors. However, IRA's use of restart- 
ing allows it to gradually improve the accuracy of a set 
of eigenstates during iteration, using a reasonably small 
number of Ritz vectors (typically a few hundred at most). 

In the extended HOSPHE code, we use the numeri- 
cal software ARPACK [2(|, which implements the IRA 
method. With this method, the number of matrix- vector 
products needed in order to reach convergence depends 
on the requested tolerance. Denoting the QRPA matrix 
A, the approximate eigenvector x and the corresponding 
approximate eigenvalue A, the iterations proceed until 
the accuracy measure \\Ax — Xx\\ [2d| is less than the re- 
quested tolerance. 

Fig. Q] shows the convergence as a function of the basis 
size when applying the method for the calculation of the 
lowest 2 + state in 214 Pb. The only truncation employed 
is the number of main oscillator shells used for the basis. 
As seen from this figure the accuracy measure is always 
lower than the requested tolerance when the iterations 
finish and the number of iterations required in order to 
reach the desired convergence increases with the size of 
the basis. 

To find the lowest eigenstates a typical choice is to start 
from a random initial guess (pivot) vector. For states 
with large transition probabilities, the number of itera- 
tions needed can however be reduced by instead starting 
from an initial pivot vector whose matrix elements are 
set to the matrix elements of the corresponding electro- 
magnetic multipole operator [l|. In the case where pair- 
ing disappears (and the numerical accuracy is high) the 
electromagnetic pivot also filters out the states that have 
an overlap with the pivot and thus removes states which 
correspond to pair addition or removal. Because of these 
advantageous features we start from an electromagnetic 
pivot in all the calculations presented. 

For the calculations presented below, a value of 17 os- 
cillator shells (N max = 16) was chosen to offer a good 
balance between accuracy and computational speed. As- 
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Figure 1: Convergence of the excitation energy and reduced 
transition probability for the lowest 2 + state in 214 Pb using the 
SLy4 interaction together with a separable Gaussian pairing 
force [lg, LL6J] . The convergence is shown as a function of the 
maximum oscillator shell N max included in the basis. The 
accuracy measure and the number of iterations needed to ob- 
tain convergence are also shown. The tolerance parameter 
which determines when to stop the Arnoldi iterations was set 
to 0.001. For N max = 16, the wall clock time was 12 min 
on a desktop workstation (Intel Core i7-2600K, 3.4GHz). In 
all cases the Krylov subspace was constructed from 100 Ritz 
vectors which was estimated to give the fastest convergence 
with N max = 16. 



suming that the calculation with 33 oscillator shells 
shown in Fig. [T] is fully converged, the truncation error 
when stopping at 17 shells amounts to 0.01 MeV for the 
energy and 0.002 (eb) for the reduced transition proba- 
bility. Using 17 shells reduces the dimension of the QRPA 
matrix to 8016 as compared to 59296 in the case of 33 
shells. With this smaller basis and using a tolerance pa- 
rameter of 0.001, the average time to calculate the lowest 
state for a nucleus in the lead isotope chain is 6.5 min 
(Intel Core i7-2600K, 3.4GHz) and the average number 
of iterations required is 2663. 



Sets of a few lowest eigenvalues can also be extracted 
and requires about the same number of iterations. For 
example, in in the case of 192 Pb, the number of itera- 
tions needed to extract 10, 20 and 30 positive energy 
eigenstates becomes 1974, 2648 and 3427 respectively. 



Interaction G„ G p 

SLy4 655 600 

SKM* 610 550 

SkX 560 530 

Table II: Strength parameters of the separable Gaussian pair- 
ing interaction in units of MeVfm' ! . For the range of the 
interaction we adopt the value a = 0.660 fm in all cases. 
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normal density-matrices are truncated during the HFB 
iterations [13]. However, in the subsequent QRPA calcu- 
lation we used non-truncated wavefunctions without any 
energy cut for the residual particle-particle interaction. 
This way of using the equivalent spectra method is there- 
fore slightly inconsistent and a better truncation recipe is 
desired. In the following we will only use the finite-range 
pairing interaction which does not need to be truncated 
and allows us to treat HFB and QRPA in a consistent 
way. 

The effect of changing the strength of the Gaussian 
pairing interaction with ±5% is shown in Fig. [5] with 
dashed lines. Both the energies and the transitions are 
sensitive to such a change, i.e. a decreese of the pairing 
lowers the excitation energies and raises the B (E2) val- 
ues. The effect is seen to be largest for N = 104 which 
is just between two magic numbers. The schematic fits 
of the pairing strengths appear to give quite reasonable 
values for the lead isotopes with an average E 2 + energy 
that agrees roughly with experiment. 



B. J w — 2 + states in Pb and Sn isotopes 



Figure 2: Excitation energies and reduced transition proba- 
bilities for Pb isotopes. Results are shown for different treat- 
ments of the pairing interaction. The dashed lines denote the 
result of changing the strength of the finite-range pairing with 
±5%. Decreasing the pairing lowers the energies and raises 
the B (E2) values. The strength parameters for the zero- 
range interaction were chosen as V n = —168 and V„ = —200 
(MeV(fm) 3 ) when the cutoff in the equivalent spectra was 
taken as 60 MeV and 20 MeV respectively. 



IV. NUMERICAL RESULTS 

A. Influence of the pairing interaction 

In order to study the influence of the pairing inter- 
action on 2 + states we compare the use of a zero-range 
delta interaction [PJ combined with a truncation in the 



equivalent spectra [l 
sian pairing force [151 . 



to the use of a separable Gaus- 
16j |. This force has a finite range 
and therefore does not need to be truncated. In order 
to obtain reasonable pairing, the pairing strengths are 
tuned to get the lowest quasiparticle energies to agree 
with the experimental gaps extracted in Ref. J21J using 
a four-point formula. The resulting parameters obtained 
for the finite-range interaction are shown in Table ITT] 

Results for the lead isotopes using the different pairing 
interactions are shown in Fig. [5] As seen in this figure 
there are fluctuations in the energies which depend on 
the choice of pairing force. Comparing the two pairing 
interactions, it appears that the finite-range interaction 
is slightly better in capturing the fluctuations of the ex- 
perimental energies. 

In the equivalent spectra method the normal and ab- 



In this work we consider three different Skyrme param- 
eterizations: SkM*, SLy4 and SkX. SKM* is based on 
the SkM parameters [23|, but has been adjusted further 
using results from fission barrier calculations |24j. The 
original SkM parameters were determined by considering 
both static ground state properties as well as some dy- 
namical properties including monopole and quadrupole 
resonances |23| . SLy4 was adjusted with special care 
taken to model neutron matter in order to facilitate the 
description of neutron rich nuclei [111 ]. The accuracy of 
QRPA based on these interactions was recently studied 
and compared to calculations based on the generator- 
coordinate method GCM [25|. It was found that SkM* 
reproduced experimental 2 + states more accurately than 
SLy4 and the QRPA results were similar to results ob- 
tained with the GCM. In addition, we also consider the 
SkX interaction which has been tuned with special fo- 
cus on reproducing single-particle states in double-magic 
nuclei [26| . 

Results for the lead isotopes using the three different 
Skyrme interactions are shown in the left hand panels of 
Fig. [31 For the double magic nucleus Pbi26, SkX gives 
the correct 2 + energy while the other two forces over- 
estimate this excitation energy. As the neutron number 
is reduced, the predictions show considerable differences. 
SLy4 gives zero energy solutions around N = 112, in- 
dicating a transition to a deformed ground state, while 
the other two interactions appear to be stiffer towards 
deformation and give more realistic results. 

The ground-state energies of Pbios and Pbn2 as a 
function of quadrupole deformation are shown in Fig. [4] 
The energy curves are calculated using delta pairing 
instead of the Gaussian pairing which means that the 
curves are slightly inconsistent with the QRPA calcula- 
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Figure 3: Excitation energies and reduced transition proba- 
bilities (B(E2;0^ s — > 2j")) for Pb and Sn isotopes. Results 
are shown for three different Skyrme parameterizations. The 
experimental values are taken from Ref. [23 |. 
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Figure 4: HFB ground-state energy as a function of 
quadrupole deformation /3 [13fl for Pbios and Pbn2. The cal- 
culation was performed using the hfbtho code [13 |. The local 
minimas are marked with symbols. Constant shifts of both 
curves have been applied in order to make them fit in the 
figure. 



tions, but the general features will be the same. As seen 
from this figure, both interactions predict spherical min- 
ima for Pbios, although the lowest minimum with SLy4 is 
the oblate one with quadrupole deformation /? = —0.19. 
The spherical minimum obtained with SkM* is stiffer 
than with SLy4 which is probably the reason why the 
2 + energy is predicted higher. As one moves to Pbn2, 
SLy4 gives the spherical point as a maximum with neigh- 
boring slightly deformed minima. In this case our QRPA 
calculation is likely to give a zero energy solution as the 
assumed spherical ground state is no longer stable with 
respect to quadrupole deformations. The 2 + energies ob- 
tained with the SkM* interaction agree rather well with 
experiment and seem to favor the prediction of stiffer 
energy surfaces. 

In the QRPA formalism, an expression for the operator 
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Figure 5: Structure of QRPA l\ states in lead isotopes using 
SkM*. 



which creates the excited states by acting on the QRPA 
ground state can be written as 



oi = £ zt w 
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In order to discuss the structure of the solutions, we 
label the kk! components of this creation operator us- 
ing the quantum numbers of the quasiparticle operators. 
As an example, if both k and k! refer to a quasiproton 
(quasineutron) in a 113/2 shell, the corresponding compo- 
nent is denoted as 7r(i 13 / 2 ) 2 (i / (*i3/2) 2 )- Indeed, in the 
limit when Z' kk , = this turns into the usual notation 
for writing two quasiprotons (quasineutrons) in the 113/2 
shell. 

The major oscillator iV osc . quantum number is not pre- 
served in the calculations, but for simplicity we will refer 
to the mixed orbitals using the harmonic oscillator order- 
ing. For example the lowest p 3 / 2 and P1/2 quasiparticle 
orbitals will be referred to as being of N osc = 1 char- 
acter, although these orbitals also contain contributions 
from higher oscillator shells. 

We denote the probability P a j, a 'j' for different compo- 
nents in the wavefunctions of the excited states by sum- 
ming contributions from the different m quantum num- 
bers as 
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Defined in this way, the largest components in the cal- 
culated 2f states in the chain of lead isotopes are shown 
in Fig. [5] As seen in this figure, the largest proton and 
neutron components in Pbi26 involve particle-hole exci- 
tations across the Z — 82 and N = 126 gaps. In the 
other isotopes, where neutron pairing is active, the 2 + 
states mainly involve neutron excitations (~ 85%) with 

the dominating component being 15-30 % v (113/2) for 
7Y = 100 — 116. The calculated neutron single-particle 
levels are shown in Figure [5] and as expected the dom- 
inating components in the 2 + states involve excitations 
among the shells close to the Fermi-level. 

It should be noted that the transition strenghts are 
calculated directly from the electromagnetic operators Q 
without any effective charges. Therefore it is the smaller 
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Figure 6: Theoretical (SkX) neutron single-particle levels. 
The left set of levels is for Ca28, the middle set for Sns2, 
and the right set is for Pbi26- Positive parity levels are shown 
with full lines and negative parity levels with dashed lines. 
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Figure 7: Same as in Fig. [3] but for N — 126 isotones. 



proton components, suppressed because of the magic pro- 
ton number, which determine the electromagnetic prop- 
erties. It is also interesting to notice that above the 126 
gap, the 2J states are composed of rather pure two-quasi- 
neutron excitations to the <7g/ 2 shell. 

For the N=126 isotones shown in Fig. 0there is a sim- 
ilar accuracy as obtained for the lead isotopes. To have a 
better idea about the amount of collectivity that should 
be present when going away from 82Pbi26, more exper- 
imental transition probabilities are clearly needed and 
experiments to measure the unknown B (E2) values for 
the isotopes 82Pbii 4! n 6 . n 8! i2o, 7sPti22,i24 and 80 Hg 12 6 
are planned [5J]. 

Results for the Sn chain are shown in the right hand 
panels of Fig. [3] As in the case of the Pb chain, the SLy4 
interaction gives some zero energy solutions, while the 
other two interactions produce dips in the excitation en- 
ergies for some isotopes, but do not reach zero. For the 
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Figure 8: Same as in Fig. [3] but for Ni and Ca isotopes. 



two double magic Sn isotopes, the excited states are cal- 
culated to be roughly even mixture of proton and neutron 
excitations, while the excited states in the semi-magic iso- 
topes mainly involve neutron excitations. With SkX, the 

largest components are v (^5/2) f° r Sns2_54, v (37/2) 
.2 



for Sn56_64 and v [h 



11/2J 



for Sng6-80- The positions of 



these shells (shown in Figure [6]) are thus important in 
order to reproduce the details of the experimental data. 



C. J w — 2 + states in Ni and Ca isotopes 

The results for the Ni chain are shown in the left hand 
panels of Fig. [8] In M28, the 2 + state is built from an 
almost equal mixture of proton and neutron excitations , 
while the 2 + states in the other nickel isotopes are dom- 
inated by neutron excitations. SkX predicts a smaller 28 
gap than the other interactions and a 2 + state in Ni 2 8 
which is lower than in experiment. SkX also predicts a 
smaller gap at 7V=32 and does not show the spike in ex- 
citation energy obtained with the other interactions for 
Ni 32 - 

QRPA calculations based on the relativistic-mean-field 
model for Ni4o [27| overpredicted the energy of the 2 + 
state by roughly three times the experimental energy. 
A suggested explanation was missing 2p-2h and higher 
order excitations among the neutrons [27| . Since the 
neutron Fermi-level is located between opposite parity 
shells, this state will be overpredicted whenever the neu- 
tron pairing goes to zero. With SLy4 we also obtain 
an overprediction, although less severe, while the other 
interactions predict excitation energies close to the ex- 
perimental value. It is interesting to note that with SkX 
the ground state is calculated to have an average gap [28| 
of A„ = 1.38 MeV and the 2 + state is obtained with both 
correct energy and transition strength. The state is built 
as a mixture of proton (~ 23%) and neutron (~ 77%) 
excitations, where the largest component is v (gQ/2) ■ 

Results for the Ca chain are shown in the right hand 
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Figure 9: Excitation energies and reduced transition proba- 
bilities (B(E3; 0^. s . — > 3J")) for the lowest 3~ states in Pb and 
Sn isotopes. Results are shown for three different Skyrme pa- 
rameterizations. The experimental data for the 3~ states are 
taken from the compilation |29l |. 



panels of Fig. [5] For Ca2o one should notice that the 
Fermi-levels for both neutrons and protons are right be- 
tween shells with opposite parity. Therefore, if pairing 
disappears, as happens with SkX, the lowest particle- 
hole excitations with positive parity are between shells 
of N osc = 2 and iV osc = 4 character, and the excita- 
tions have to bridge an energy gap of around 15 MeV, 
which pushes the predicted 2 + state up to a high energy. 
With SkM* (SLy4) there is some pairing remaining and 
a low-lying 2 + state can be constructed with a dominat- 
ing 7r (0^3/2) (tt (77/2) ) quasiparticle component. This 
state is likely to have an average particle number which 
is somewhat wrong, but in this work we only remove the 
excitations being in the wrong nucleus if pairing vanishes 
completely. Notice also that since the transition opera- 
tor is of particle-hole type, excitations corresponding to 
addition or removal of two particles gives zero for the 
transition strength. 

In general, the 2 + states in the Ca chain are predicted 
to have too little collectivity compared to experiments. 
The states in Ca22-26 are predicted to be rather pure 
excitations within the /V/2 shell. For exanrple, with the 

SkX interaction the component of v (/7/2J is 95, 90, 
and 82% for the excitations in Ca22-26 respectively. In 
order to induce more collectivity, it appears likely that 
proton two-particle-two-hole excitations across the Z = 
20 gap must be explicitly included which goes beyond the 
present QRPA treatment. 



trons are added, negative parity states can be made by 
exciting particles between the neighboring positive parity 
gg/2 shell and negative parity j'15/2 shell, located above 
the 126 gap (see Figure [5]). Therefore the 3~ state in 
Pt>r28 can be built at a low cost mainly from neutron 
excitations, which explains the dip in the experimental 
energy. As seen in Fig. |H1 SkX reproduces the experimen- 
tal lowest 3~ energies and B(E3) values almost perfectly 
except for the dip in the transition probability seen for 
Pbi28- The forces SkM* and SLy4 have a less perfect 
overall agreement, but SkM* agrees with experiment in 
the case of Pbi28- 

In Pbioo with SkX, the excitation operator for the low- 
est 3~ state consists of 17 % proton excitations, where 
the largest component is only 6 % ir (d.3/2 )* (^9/2 )* ail d the 
main neutron component is 62 % v(fr/2) 1 (H3/2) 1 ■ The 
3~ states of the neutron deficient Pb isotopes have many 
small proton components contributing. When the neu- 
tron number increases, proton excitations become more 
donrinant, but with fewer components contributing. At 
#=124, proton excitations constitute 48 % with one 
donrinant component, 29 % 7r(d 3 /2) 1 (/ig/2) 1 j while the 
largest neutron component is 13 % v{jpz/2) 1 (<?9 /2) 1 ■ Thus, 
with QRPA, we see a transition from strong proton con- 
figuration mixing in Pbioo to strong neutron configura- 
tion mixing in Pbi24. 

The lowest experimental and calculated 3~ states of 
tin isotopes are shown in the right hand panels of Fig. [9] 
The energies are reproduced almost perfectly in the re- 
gion iV=60-72, using SkX and SkM* forces. Closer to 
the 7Y=82 shell closure, the predicted energies become 
too high. The experimental B(E3) values are repro- 
duced roughly by all three Skyrme forces, but none of 
the three forces gives a truly accurate description of 
the finer details. The excitations are dominated by the 
i/ ( ( ^5/2) 1 (^ir/2) 1 configuration which has a component of 
70 % in Sn 52 , and gradually goes down to 50 % in Sn 8 o- 
The second largest excitation is 3% ^(<?9/2) 1 (^n/2) 1 m 
Sn52, but starting from Sn54 the second largest excita- 
tion is i / (37/2) 1 (^n/2) 1 an< i its anrplitude grows steadily 
from 5% to 10% in Snso. Other neutron excitations be- 
tween the vh\\/2 subshell and the N osc = 4 orbitals are 
excluded because of angular momentum selection rules. 
The proton fraction of the transitions is an almost con- 
stant 10 % and is composed of many small-amplitude 
excitations. The lowest energies and the largest transi- 
tions strengths are obtained for 7Y=66 when the neutron 
Fermi-level is close to the vhn/2 intruder shell and neg- 
ative parity excitations correspond to a low energy cost. 



D. J n = 3 states in Pb and Sn isotopes 

The calculated 3~ states for lead isotopes are shown in 
the left hand panels of Fig. [5] The most striking feature 
of the data is the dip in excitation energy seen when go- 
ing from Pbi26 to Pb^s- In Pt>i26 the 3~ state is created 
by an roughly equal amount of neutron and proton exci- 
tations across the 82 and 126 gaps. When two more neu- 



E. J 71 — 3 states in Ni and Ca isotopes 

For the Ni isotopes shown in the left hand panels of 
Fig. [10] we also show experimental data for the second 
3~ states in Ni28-32 taken from [30J]. For some reason 
the calculations for Ni28-32 agree better with the sec- 
ond 3~ states. Especially the lowest 3~ state in N128 
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Figure 10: Same as in Fig. [9l but for Ni and Ca isotopes. 
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Figure 11: Leading two-quasiparticle components of QRPA 
3j~ states in Ca isotopes using SkX. 



the largest components in Ca4Q. 



is predicted about 3 MeV too high in energy. It should 
be noted that this state has not been clearly identified 
as 3~ in experiments [29j, contrary to the 3 _ states in 
the other Ni isotopes. One should also note that with 
Skyrme type interactions one usually neglects proton- 
neutron pairing and part of the isovector particle-hole 
interaction which may have an influence around the N=Z 
line. With SkX, the lowest theoretical 3~ state in M28 is 
calculated to be of isoscalar type. This state is built from 
28 % v(d 3/2 ) 1 {p 3/2 ) 1 and 22 % 7r(d 3/2 ) 1 (p3 /2 ) 1 along with 
smaller probability excitations to the orbitals of N osc = 4 
character. 

For the semi-magic nickel isotopes with 7V=30-48, 
the 3 _ states are dominated by neutron excitations, 
l/ (P3/2) 1 (99/2) 1 which decrease from 75 % in iV=30 to 
16 % in 7Y=48 and v(f5/2) 1 (g9/2) 1 > which increase from 
2 % in A=30 to 74 % in JV=48. The proton components 
are small since it is more favorable to excite neutrons 
when the proton Fermi- level is just above /V/2- The dou- 
ble magic M50 has a different structure than the semi- 
magic isotopes, the dominating component being 36 % 

v(Pi/2) 1 (d 5 /2) 1 - 

The results for calcium isotopes are shown in the right 
hand panels of Fig. [TU] In this case, the lowest exper- 
imental 3~ energies are best reproduced by SkM* and 
SkX while SLy4 gives too high energies. However, in 
general the finer details of the 3J" energies between N=20 
and iV=26 are not reproduced. Especially for the N = Z 
nucleus Ca2o , the energy calculated with SkM* and SkX 
becomes too low. 

The leading wave function components of the calcium 
isotopes are shown in Fig. [TTJ As seen in this figure, the 
3~ state in 2oCa2o is composed of a fairly even mixture 
of proton and neutron excitations from the d 3 « and Si / 2 
shells below the 20 gap to the / 7 / 2 shell just above the 
gap (see Figure [5]). Going towards Ca28 the neutron exci- 
tations become suppressed as the Fermi-level reaches the 
middle of the N osc = 3 shell and proton excitations start 
to dominate. When more neutrons are added, neutron 
excitations to the g 9 / 2 shell start to appear and become 



V. SUMMARY AND CONCLUSIONS 

An iterative method for the solution of the QRPA 
equations which avoids the construction of the large 
QRPA matrix was employed for the calculation of low- 
lying vibrational states. The method uses the Implicitly 
Restarted Arnoldi approach for the solution of the non- 
hermitian eigenvalue problem. In this approach, only 
the action of the matrix on a Ritz vector is needed. As 
demonstrated, this can be expressed in terms of the fields 
generated by the transitional densities corresponding to 
the Ritz vector. Our study shows that the method is nu- 
merically stable and typically requires a few thousand it- 
erations in order to produce well converged lowest eigen- 
states. 

The new solution method was applied to the calcula- 
tion of excitation energies and decay rates of the first 2 + 
and 3~ vibrational states in a set of spherical even-even 
nuclei. The calculations were performed using three dif- 
ferent Skyrme interactions together with a finite-range 
pairing force. Overall a quite reasonable agreement with 
experimental data was obtained. The main difficulties 
seem to be in the description of 2 + states right between 
two magic neutron numbers where the different interac- 
tions tend to give different results and even zero energy 
solutions. 

Difficulties were also observed for the Ca isotopes 
where all the interactions gave too little collectivity. 
These difficulties are probably related to the limitations 
of the QRPA method itself, as it only includes two- 
quasiparticle excitations. However, since our method is 
rather computationally inexpensive, it may become prac- 
tical to consider extensions of QRPA, for example higher- 
order QRPA approaches or boson expansion methods, 
which allow to treat more complicated excitations that 
could improve the results. 

Because of the low numerical cost and low memory 
requirements, the method appears promising for appli- 
cations to deformed nuclei where the dimensions beconre 
substantially larger. Indeed, the methods of this work 



10 



are quite analogous to the ones used in the nuclear Shell 
Model community, where quite similar iterative methods 
have been used for large-dimensional hermitean eigen- 
value problems. However, contrary to the Shell Model 
where the dimensions increase exponentially and multi- 
major shell calculations for heavy nuclei are almost im- 
possible, the QRPA method stays tractable. Therefore, 
as a next step, we will implement the method in a code 
able to treat nuclei with deformed ground states. 

For spherical nuclei, the speed of the iterative method 
opens the possibility to directly compute low-lying states 
and include them as part of the observables used when 
fitting the parameters of new improved Skyrme interac- 
tions. However, care must be taken to analyse the struc- 
ture of the included states in order to ensure that a QRPA 



description is compatible with the experimental states. 
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